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Abstract. Euler-Lagrange equations and variational integrators are developed for Lagrangian mechanical systems evolving 
on a product of two-spheres. The geometric structure of a product of two-spheres is carefully considered in order to obtain global 
equations of motion. Both continuous equations of motion and variational integrators completely avoid the singularities and 
complexities introduced by local parameterizations or explicit constraints. We derive global expressions for the Euler-Lagrange 
equations on two-spheres which are more compact than existing equations written in terms of angles. Since the variational 
integrators are derived from Hamilton's principle, they preserve the geometric features of the dynamics such as symplecticity, 
momentum maps, or total energy, as well as the structure of the configuration manifold. Computational properties of the 
variational integrators are illustrated for several mechanical systems. 
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1. Introduction. The two-sphere is the set of aU points in the Euchdean space M.^ which are a 
unit distance from the origin. It is a two dimensional manifold that is locally diffeomorphic to R^. Many 
classical and interesting mechanical systems, such as a spherical pendulum, a double spherical pendulum, 
and magnetic models, evolve on the two-sphere or on a product of two-spheres. In this paper, we derive 
Euler-Lagrange equations on configuration spaces of the form (S^)", for a positive integer n. We also develop 
geometric numerical integrators referred to as discrete Euler-Lagrange equations or variational integrators 
on (§2)". 

In most of the literature that treats dynamic systems on (§^)", either 2n angles or n explicit equality 
constraints enforcing unit length are used to describe the configuration of the system [21 [15] . These descrip- 
tions involve complicated trigonometric expressions and introduce additional complexity in analysis and 
computations. In this paper, we focus on developing continuous equations of motion and discrete equations 
of motion directly on (§^)", without need of local parameterizations, constraints, or reprojections. This 
provides a remarkably compact form of the equations of motion. 

Geometric numerical integrators are numerical integration algorithms that preserve the geometric struc- 
ture of the continuous dynamics, such as invariants, symplecticity, and the configuration manifold [5]. Con- 
ventional numerical integrators construct a discrete approximation to the flow using only information about 
the vector field, and ignore the physical laws and the geometric properties inherent in the differential equa- 
tions 17J. Consequently, they do not preserve important characteristics of the dynamics of the continuous 
equations of motion. In contrast, variational integrators are constructed by discretizing Hamilton's principle, 
rather than discretizing the continuous Euler-Lagrange equation [18l [16] . Since they are developed by using 
a discrete version of a physical principle, the resulting integrators have the desirable property that they are 
symplectic and momentum preserving, and they exhibit good energy behavior for exponentially long times. 

Geometric numerical integration on has been studied in [121 [HI [13] • The two-sphere is a homogeneous 
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manifold; the special orthogonal group SO (3) acts transitively on S^, and Lie group methods [7] can be 
adapted to generate numerical flows on S^. 

In this paper, we study Lagrangian mechanical systems on (S'^)". Thus, it is desirable to preserve the 
geometric properties of the dynamics, such as momentum map, symplecticity, and total energy, in addition 
to the structure of the configuration manifold llj. We combine the approaches of geometric integrators on 
homogeneous manifolds and variational integrators to obtain variational integrators on (S^)" that preserve 
the geometric properties of the dynamics as well as the homogeneous structure of the configuration manifold 
(§^)" concurrently. 

The contributions of this paper can be summarized in two aspects. In the continuous time setting, the 
global Euler-Lagrange equations on (§^)" are developed in a compact form without local parameterization 
or constraints. This provides insight into the global dynamics on (§^)", which is desirable for theoretical 
studies. As a geometric numerical integrator, the discrete Euler-Lagrange equations on (§^)" are unique in the 
sense that they conserve both the geometric properties of the dynamics and the manifold structure of (S^)" 
simultaneously. The exact geometric properties of the discrete flow not only generate improved qualitative 
behavior, but they also provide accurate and reliable computational results in long-time simulation. 

This paper is organized as follows. Lagrangian mechanics on is described in Section [51 Variational 

integrators on (S^)" are developed in Section [31 Computational properties are illustrated for several me- 
chanical systems, namely a double spherical pendulum, an n-body problem on a sphere, an interconnected 
system of spherical pendula, pure bending of an elastic rod, a spatial array of magnetic dipoles, and molecular 
dynamics that evolves on a sphere. 

2. Lagrangian mechanics on (§^)". In this section, continuous equations of motion for a mechanical 
system defined on (§^)" are developed in the context of Lagrangian mechanics. It is common in the published 
literature that the equations of motion are developed by using either two angles or a unit length constraint 
to characterize Any description with two angles has singularities, and any trajectory near a singularity 
experiences numerical ill-conditioning. The unit length constraint leads to additional complexity in numerical 
computations. We develop global continuous equations of motion without resorting to local paramctcrizations 
or constraints. To achieve this, it is critical to understand the global characteristics of a mechanical system 
on (S^)". This section provides a good background for understanding the theory of discrete Lagrangian 
mechanics on (S^)" to be introduced in the next section. 

The two-sphere is the set of points that have the unit length from the origin of M'^, i.e. §^ = {(76 
M'^ I g • (7 = 1}. The tangent space TgS^ for (7 G §^ is a plane tangent to the two-sphere at the point q. Thus, 
a curve 5 : M ^ §^ and its time derivative satisfy q ■ q = Q. The time-derivative of a curve can be written as 



where the angular velocity (jj G K."^ is constrained to be orthogonal to q, i.e. q ■ uj — Q. The time derivative 
of the angular velocity is also orthogonal to g, i.e. q ■ to — Q. 

2.1. Euler-Lagrange equations on (§^)". We consider a mechanical system on the configuration 
manifold §^ x • • • x §^ = (§^)". We assume that the Lagrangian L : r(§^)" E is given by the difference 
between a quadratic kinetic energy and a configuration-dependent potential energy as follows. 



q ^ uj X q, 



(2.1) 




(2.2) 



Lagrangian Mechanics and Variational Integrators on Two-Spheres 



3 



where {qi,qi) G TS^ for i G {1, . . . ,n}, and Af^ G R is the i,j-th element of a symmetric positive definite 
inertia matrix M G M"^" for i,j G {I,.. .,n}. The configm-ation dependent potential is denoted by V : 

The action integral is defined as the time integral of the Lagrangian, and the variation of the action 
integral leads to continuous equations of motion by applying Hamilton's principle. These are standard 
procedures to derive the Euler-Lagrange equations. The expression for the infinitesimal variation of qi G §^ 
should be carefully developed, since the configuration manifold is not a linear vector space. As in (|2.ip . the 
infinitesimal variation of qi can be written as a vector cross product. 



Sqi = X gj. 



(2.3) 



where £,i G R'^ is constrained to be orthogonal to qi, i.e. • = 0. From this, the expression for the 
infinitesimal variation of q,; is given by 



Sqi = 6 X gi + 6 X qi. 

These expressions are the key elements to obtaining global equations of motion on (§^)". 
The variation of the Lagrangian can be written as 



(2.4) 



SL= ^ Sqi ■ Mijqj - ^ Sqi 



dV 



i=l 



where the symmetric property Mij ~ Mji is used. Substituting (|2.3p and (|2.4p into this, and using the vector 
identity [a x h) ■ c — a ■ {b y. c) for any a, 5, c G M^, we obtain 

SL= ^ i,- (g, X M.^qj) + ■ {q, x A'hjqj) - ■ Iqi >^ g-j ■ 

Let 25 be the action integral defined as 25 = L{qi, . . . , g„, gi, . . . , g„) dt. Using the above equation and 
integrating by parts, the variation of the action integral is given by 



S&= ^ • (9i X AhjCij + g, X Mijiij] 



T n „T 

.=1 -'o 



{qi X 2_^M,jqj) + qi x — 
3=1 '^'^^ 



From Hamilton's principle, (525 — for any ^i vanishing at t = 0,T. Since S,i is orthogonal to qi, the 
continuous equations of motion satisfy 



, v-^ , ^ ... dV , , 

X 2^ Afygj) + gj X — = Ci{t)qi 

3 = 1 "^^^ 



(2.5) 



for some scalar valued functions Ci(i) for i G {1, n}. Taking the dot product of (|2.5p and qi implies that 
= Ci(t)||gi||^ = Ci(t), which is to say that the scalar valued functions are uniformly zero. Now we find an 
expression for qi. Since the left hand side expression is perpendicular to qi, it is zero if and only if its cross 
product with qi is zero. Thus, we obtain 



X (g, X ^ M^jiij) + g, X ( g, x 

3 = 1 



dV 



0. 



(2.6) 
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From the vector identity a x (6 x c) = (a • c)b — (a • b)c for any a,b,c G R^, we have 



Qi X {qi X iji) = {qi ■ qi)qi ~ {qi ■ qi)qi, 

= -{4i ■ qi)qi - 

where we use the properties jiiqi ■ qi) — qi ■ qi + qi ■ qi — and qi ■ qi — 1. Substituting these into (|2.6p . we 
obtain an expression for qi, which is summarized as foUows. 

Proposition 1 Consider a mechanical system on (S^)" whose Lagrangian is expressed as S2.^) . The con- 
tinuous equations of motion are given by 



Maqi = qiX {qix'^ M^jqj) - {qi ■ qi)Miiqi + q^ x [qi x — j 



(2.7) 



for z G {1, . . . , 7i}. Equivalently, this can be written in matrix form as 



^11^3x3 


-Muqiqi ■ ■ 


• -Minqiqi 




qi 






M22I3X3 


■ -M2nq2q2 




h 




-MniqnCjn 


-Mn2qnqn ' ' 






qn_ 





-{qi-qi)Miiqi+qf§^ 
-{q2-q2)M22q2 + ql^ 



-{qn ■ qn)Mn„q„ + ql§^ 



(2.8) 



where the hat map ' : M.^ ^ MP^^ is defined such that ab — a x b for any a, 6 G M'^. 

Since qi = uji x qi for the angular velocity uii satisfying qi ■ uji = 0, we have 

(ji ^ uJi X qi + X {uJi X qi), 
^ uJi X qi - (tjj • uJt)qt. 

Substituting this into (|2.5p and using the fact that qi ■ Coi — 0, we obtain continuous equations of motion in 
terms of the angular velocity. 

Corollary 1 The continuous equations of motion given by |i?.7| j can be written in terms of the angular 
velocity as 



dV 



MiiLOi = ^ {Mijqi X {qj X ujj) + Mij{ujj ■ UJj)qi X qj) - <?,; X — 

cfqi 



qi =Ui X qi 

for i G {1, . . . , n}. Equivalently, this can be written in matrix form as 



(2.9) 
(2.10) 



Mii/3x3 


-Mi2qiq2 ■ ■ 


• -Minqiqn 








-M2ig29i 


M22I3X3 


■ -M2nq2qn 




LO2 




-M„iq„qi 


-Mn2qnq2 ' ' 











. E"=l ^lnj{^j ■ UJj)qnqj - qn§^ . 



(2.11) 



Lagrangian Mechanics and Variational Integrators on Two-Spheres 



5 



Equations (|2.7p - (|2.1ip are global continuous equations of motion for a mechanical system on (§^)". They 
avoid singularities completely, and they preserve the structure of T(§^)" automatically, if an initial condition 
is chosen properly. These equations are useful for understanding global characteristics of the dynamics. In 
addition, these expressions are dramatically more compact than the equations of motion written in terms of 
any local parameterization. 

We need to check that the 3n x 3n matrices given by the first terms of (|2.8p and (|2.1ip are nonsingular. 
This is a property of the mechanical system itself, rather than a consequence of the particular form of 
equations of motion. For example, when n = 2, it can be shown that 

M11/3X3 -Mi2qiq2 
-Mi2q2qi M22/3X3 

= Mi\M|2(MiiM22 - Mf2(9i ■ q2 f){MiiM22 - Affa). 

Since the inertia matrix is symmetric positive definite, Mn, M22 > 0, M11M22 > M12, and from the Cauchy- 
Schwarz inequality, {qi ■ (72)^ £ (91 • qi){q2 ■ 92) — 1- Thus, the above matrices are non-singular. One may 
show a similar property for n > 2. Throughtout this paper, it is assumed that the 3n x 3n matrices given 
at the first terms of (12. 8|) and (|2.11[) are nonsingular. Under this assumption, the Legendre transformation 
given in the next subsection is a diffeomorphism; the Lagrangian is hyperregular. 

2.2. Legendre transformation. The Legendre transformation of the Lagrangian gives an equivalent 
Hamiltonian form of equations of motion in terms of conjugate momenta if the Lagrangian is hyperregular. 
Here, we find expressions for the conjugate momenta, which are used in the following section for the discrete 
equations of motion. For q^ G EP, the corresponding conjugate momentum pi lies in the dual space T*.S^. 
We identify the tangent space Tg.S^ and its dual space T*.E>'^ by using the usual dot product in R^. The 
Legendre transformation is given by 

Pi ■ 6qi = D^,L(9i, . . . ,g„,gi, . . . ,g„) • 

n 

which is satisfied for any Sqi perpendicular to qi. Here D^.L denotes the derivative of the Lagrangian with 
respect to qi. The momentum pi is an element of the dual space identified with the tangent space, and the 
component parallel to qi has no effect since Sqi ■ qi — 0. As such, the vector representing pi is perpendicular 
to qi, and Pi is equal to the projection of X]j=i ^'^ij^j onto the orthogonal complement to qi, 

n n 

= ^i^Iij'ij " ill ■ Mijqj)qi) = ^{{qi ■ qi)M^jqj - {qi ■ ]Vh,q,)qi), 

n 

= Mi.qi - gi X (gi X ^ A%«i)- (2-12) 

3. Variational integrators on (§^)". The dynamics of Lagrangian and Hamiltonian systems on (§^)" 
have unique geometric properties; the Hamiltonian flow is symplectic, the total energy is conserved in the 
absence of non-conservative forces, and the momentum map associated with a symmetry of the system is 
preserved. The configuration space is a homogeneous manifold. These geometric features determine the 
qualitative dynamics of the system, and serve as a basis for theoretical study. 



det 



Mii/sxs -Mi2qiqi 
-Mi292g2 M22/3X3 



det 
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Conventional numerical integrators construct a discrete approximation of the flow using only information 
about the vector field. Other than the direction specified by the vector field, they completely ignore the 
physical laws and the geometric properties inherent in the differential equations |17| . For example, if we 
integrate (|2.1ip by using an explicit Runge-Kutta method, the unit length of the vector g^, and the total 
energy are not preserved numerically; we will see this later in this paper. 

Numerical integration methods that preserve the simplecticity of a Hamiltonian system have been stud- 
ied [20]. Coefficients of a Runge-Kutta method can be carefully chosen to satisiy a simplecticity criterion 
and order conditions to obtain a symplectic Runge-Kutta method. However, it can be difficult to construct 
such integrators, and it is not guaranteed that other invariants of the system, such as the momentum map, 
are preserved. Alternatively, variational integrators are constructed by discretizing Hamilton's principle, 
rather than by discretizing the continuous Euler-Lagrange equation [18[ 116] . The key feature of variational 
integrators is that they are derived by a discrete version of a physical principle, so the resulting integra- 
tors satisfy the physical properties automatically in a discrete sense; they are symplectic and momentum 
preserving, and they exhibit good energy behavior for exponentially long times. Lie group methods are 
numerical integrators that preserve the Lie group structure of the configuration space [7]. Recently, these 
two approaches have been unified to obtain Lie group variational integrators that preserve the geometric 
properties of the dynamics as well as the Lie group structure of the configuration manifold [10] . 

The two-sphere is a homogeneous manifold. It does not have a Lie group structure by itself, but instead, 
the special orthogonal group, S0(3) = {F £ R^^"^ \F'^F = /3x3,detF = 1}, acts on in a transitive 
way; for any 91,(72 G S^, there exist F e SO (3) such that (72 ~ Fqi. If a group acts transitively on a 
manifold, a curve on the manifold can be represented as the action of a curve in the Lie group on an initial 
point on the manifold. As such. Lie group methods can be applied to obtain numerical integration schemes 
for homogeneous manifolds [T51 [T31 [H] . However, it is not guaranteed that these methods preserve the 
geometric properties of the dynamics. In this paper, we focus on a Lagrangian mechanical system evolving 
on the homogeneous manifold, (§^)" by extending the method of Lie group variational integrators pT| llOj. 
The resulting integrator preserves the dynamic characteristics and the homogeneous manifold structure 
concurrently. 

3.1. Discrete Euler-Lagrange equations on (§-^)". The procedure to derive discrete Euler-Lagrange 
equations follows the development of the continuous time case; the tangent bundle is replaced by a cartesian 
product of the configuration manifold, a discrete Lagrangian is chosen to approximate the integral of the 
Lagrangian over a discrete time step, and the variation of the corresponding discrete action sum provides 
discrete Euler-Lagrange equations, referred to as a variational integrator. The discrete version of the Lcgendre 
transformation yields the discrete equations in Hamiltonian form. 

Let the number of timesteps be A'^, with constant timesteps h > 0. A variable with subscript k denotes 
the value of variable a,tt — kh. Define a discrete Lagrangian Ld : (S'^)" x (S^)" M such that it approximates 
the integral of the Lagrangian given by (|2.2|) over a discrete time step 

1 " h h 

- •,gn,,gi,+i, - •■ ,9n,+i) = ^ XI -*J • fefe+i - g^J - i^Vk - -Vk+i, (3.1) 

where Vk denotes the value of the potential at the k-ih step, i.e. T4 = V{qi^^, . . . ,qnk)- As given in (|2.3|) . 
the infinitesimal variation of qi^ is written as 



(3.2) 
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where G R"^ is constrained to be orthogonal to g,,., i.e. ■ qi^, = 0. The variation of the discrete 
Lagrangian can be written as 

Substituting p.2p into p.3p . and using the vector identity {a x b) ■ c — a ■ {b x c) for any a, 6, c £ M'^, we 
obtain 

1 " 



Let be the discrete action sum defined as 0d — '^k=o ^dk i which approximates the action integral as 
the discrete Lagrangian approximates a piece of the action integral over a discrete time step. The variation 
of the action sum is obtained by using (13. 4|) . Using the fact that ^i^. vanish at fc = and k — N, we can 
reindex the summation, which is the discrete analog of integration by parts, to yield 



AT-l .„ 



fc=i 1=1 



1 " 



X ^ 



From discrete Hamilton's principle S&d = for any perpendicular to qi^ . Using the same argument given 
in (|2.5p . the discrete equations of motion are given by 

1 -A dVk 

^(9*. X 2^ M^j{-qj,^, + 2q.j^ - J) - hq^^ X — = (3.5) 

for i e {1, . . . n}. In addition, we require that the unit length of the vector qi^ is preserved. This is achieved 
by viewing §^ as a homogeneous manifold. Since the special orthogonal group SO (3) acts on transitively, 
we can define a discrete update map for qi^ as 

for Fi^ G SO (3). Then, the unit length of the vector qi is preserved through the discrete equations of motion, 
since qik+i ■ qik+i — iTk-^Z-^iklik = 1- These results are summarized as follows. 

Proposition 2 Consider a mechanical system on (§^)" whose Lagrangian is expressed as i2.2\) . The discrete 
equations of motion are given by 

Muqik X Fi^q,^ + g^, x ^ Mij{Fj, ~ hx3.)qjk = ft. x ^ M,j{qj^ - gj,_ J - h^qt^ x = 0, (3.6) 

'l^k + l ^ F^k1^k (3-7) 

for z £ {1, . . . n}. For given {qik_i,qik)> we solve iS. b]) to obtain Fi^ G S0(3). Then, qi^^-^ is computed by 
Jg. 7[ ). This yields a discrete flow map {qik_i,qik) (ftkiftt+i); '^'^^ t^^^ process is repeated. 
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3.2. Discrete Legendre transformation. We find discrete equations of motion in terms of the an- 
gular velocity. The discrete Legendre transformation is given as follows |16j . 



1 " 



hdVk 



which can be directly obtained from p.3p . This is satisfied for any 6qii_ perpendicular to qi,, . Using the same 
argument used to derive (j2.12p . the conjugate momenta pi^. is the projection of the expression in brackets 
onto the orthogonal complement of qt^. . Thus, we obtain 



1 " h 



dVk 



Comparing this to (j2.12p , substituting qt^, ~ lOi^, x qi^ , and rearranging, we obtain 



n \ ^ h 



dVk 

dq.,^ 



0. 



Since the expression in the brackets is orthogonal to g^^. , the left hand side is equal to zero if and only if the 
expression in the brackets is zero. Thus, 



" 1 " h 



dVk 



(3.8) 



This provides a relationship between {qi^^uji^) and {Qi^.^ Qik+i)- Comparing this with (|3.5[) . we obtain 



(3.9) 



which provides a relationship between (g^^ , w^j. ) and {qi^_-^ j "Zu )■ Equations p.Sp and p.9p give a discrete flow 
map in terms of the angular velocity; for a given (g^j. , w^j. ) , we find {qikilik+i) by using p.Sp . Substituting 
this into (|3.9p expressed at the fc + 1th step, we obtain (g^^^^ , w^^^^ ). This procedure is summarized as 
follows. 



Corollary 2 The discrete equations of motion given by k3. 6\) and {3. 7| ) can he written in terms of the angular 
velocity as 



M^^qik X Fi^qi^ + q,^ x ^ Afjj (i^j, - I3x3)qjk = M,,huJi^ - {qi^ x ^ Afy(gj, x /iw^J) - — g,, 



(3.10) 
(3.11) 

Afi.w,,^, - (gu-+i X ^My(gj,^^ X Wj.^J) = t(9»;=+i x ^ A/y(gj,^^ - g^J) - -g,;,^^ x (3.12) 



1 

h' 
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for z G {1, . . . , n}. Equivalently, i3.1S\) can be written in a matrix form as 



-M2ig2fc+i<?i,+i M22/3X3 



-M2«g2fc+i'7nfc + i 



(3.13) 



For a given [qi^^uji^), we solve i3.10\) to obtain Fi^ G S0(3). Then, Qi^^-^ and uji^^-^ are computed by 
113.11]) and i3.13\) . respectively. This yields a discrete flow map in terms of the angular velocity (qi^. , ) i-^ 
((jj^^j , Wij.^^ ), and this process is repeated. 



3.3. Computational approach. For the discrete equations of motion, we need to solve (|3.6p and 
p.lOp to obtain _Fi^ G SO (3). Here we present a computational approach. The implicit equations given by 
and p.lOp have the following structure. 



(3.14) 



for i G {1, . . . ,n}, where My G M, qi G G M"^ are known, and we need to find Fi G SO (3). We derive an 

equivalent equation in terms of local coordinates for Fi . This is reasonable since Fi represents the relative 
update between two integration steps. Using the Cayley transformation [3T], Fi G S0(3) can be expressed 
in terms of fi G M'^ as 

Fi — {hx3 + fi){hx3 — fi) ^, 



1 



{{I- f^■ f^)hx^+2f,fI + 2f,). 



The operation Fiqi can be considered as a rotation of the vector qi about the direction fi with rotation angle 
2tan^^ ll/i||- Since the rotation of the vector qi about the direction qi has no effect, we can assume that fi 
is orthogonal to g^, i.e. fi ■ qi = 0. Under this assumption, Fiqi is given by 



Fiqi = 



1 



1 + / . • /, 



((1 - /. • fi)q^ + 2/,g,)- 



(3.15) 



Thus, we obtain 



2 2 
qi X Fiq^ = — — — -q^ X (/, X qi) = —— — - 

^ y Ji ' Ji i- t Ji ■ Ji 

2 

{Fj - /3x3)gj = ^TVjPTj^'^^^'^ ^ 
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where we use the property, qifi - 

2Mii/3x3 
2Af21g2(gl+gl/f ) 



2M„ig„(gi+gi/i^) 
l + /l7l 



Qi ^ fi — ~fiQi- Substituting these into (j3.14p . we obtain 



2A/l2<?l(<?2+g2/2^) 
1 + /2-/2 
2Af22J3x3 
1 + /2-/2 



2M„29„(<?2+g2/2^) 
1 + /2-/2 



2A/i„gi(g„+9„/J) 

2A-/2„g2(g„+g„/^') 

!+/„■/„ 



2Af„„/3x 





7r 








/2 








_/n_ 




_dn_ 



(3.16) 



which is an equation equivalent to (j3.14p . written in terms of local coordinates for Fi using the Cayley 
transformation. Any numerical method to solve nonlinear equations can be applied to find fi. Then, Fiqi is 
computed by using ()3.15|) . In particular, ()3.16|) is written in a form that can be readily applied to a fixed 
point iteration method ^SJ- 

If there are no coupling terms in the kinetic energy, we can obtain an explicit solution of ()3.14p . When 



Mij = for i 7^ j, (|3.16p reduces to 



1 + /. • / 



Using the identity, i^l'^^" g = sin 20 for any G M, it can be shown that the solution of this equation is given 
by fi = tan (i sin^^(||di|| /Ma)) -[t^. Substituting this into (|3.15p and rearranging, we obtain 





d., 

















1/2 



Using this expression, we can rewrite the discrete equations of motion given in (j3.10p - (|3.13p in an explicit 
form. 



Corollary 3 Consider a mechanical system on (§^)" whose Lagrangian is expressed as i2. S\) where Mij = 
for i ^ j, i.e. the dynamics are coupled only though the potential energy. The explicit discrete equations of 
motion are given by 



huji 



2M,:, 



X g., + 1 



dVk 



h 



dVk 



2Af„ 



1/2 



k+1 



2M,, 



^q^, 2Mi, 



(3.17) 
(3.18) 



for i £ {1, 



.,n 



}• 



3.4. Properties of variational integrators on (§^)". Since variational integrators are derived from 
the discrete Hamilton's principle, they are symplectic, and momentum preserving. The discrete action sum 
can be considered as a zero-form on (§^)" x (§^)" which maps the initial condition of a discrete flow satisfying 
the discrete Euler-Lagrange equation to the action sum for that trajectory. The simplecticity of the discrete 
flow follows from the fact the the iterated exterior derivative of any differential form is zero. If the discrete 
Lagrangian exhibits a symmetry, the corresponding momentum map is preserved since by symmetry, the 
variation of the discrete Lagrangian in the symmetry direction is zero, which in combination with the discrete 
Euler-Lagrange equations, implies a discrete version of Noether's theorem. Detailed proofs for the symplectic 
property and the momentum preserving property can be found in [16\. The total energy oscillates around 
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its initial value with small bounds on a comparatively short timescale, but there is no tendency for the mean 
of the oscillation in the total energy to drift (increase or decrease) over exponentially long times [4] . 

The variational integrators presented in this paper preserve the structure of (§^)" without need of 
local parameterizations, explicit constraints or reprojection. Using the characteristics of the homogeneous 
manifold, the discrete update map is represented by a group action of SO (3), and a proper subspace is 
searched to obtain a compact, possibly explicit, form for the numerical integrator. As a result, the following 
numerical problems are avoided: (i) local parameterizations yield singularities; (ii) numerical trajectories 
in the vicinity of a singularity experiences numerical ill-conditioning; (iii) unit length constraints lead to 
additional computational complexity; (iv) reprojection corrupts the numerical accuracy of trajectories [51ll3j. 

It can be shown that these variational integrators have second-order accuracy as the discrete action sum 
is a second-order approximation of the action integral. Higher-order integrators can be easily constructed 
by applying a symmetric composition method [23j . 

3.5. Numerical examples. The computational properties of variational integrators on and 
explicit Runge-Kutta methods are compared for several mechanical systems taken from variety of scientific 
areas, namely a double spherical pendulum, an n-body problem on a sphere, an interconnected system of 
spherical pendula, pure bending of an elastic rod, a spatial array of magnetic dipoles, and molecular dynamics 
that evolves on a sphere. 

Example 1 (Double Spherical Pendulum) A double spherical pendulum is defined by two mass par- 
ticles serially connected to frictionless two degree-of-freedom pivots by rigid massless links acting under a 
uniform gravitational potential. The dynamics of a double spherical pendulum has been studied in 15], and 
a variational integrator is developed in [22] by explicitly using unit length constraints. 

Let the mass and the length of the pendulum be mi, m2, h, h £ respectively, and let 63 = [0, 0, 1] G 
be the direction of gravity. The vector qi G represents the direction from the pivot to the first mass, 
and the vector q2 G represents the direction from the first mass to the second mass. The inertia matrix 
is given by Mn — {mi + m2)ll, M12 — TO2^i^2, and M22 — 7712^2- The gravitational potential is written as 
92) — —{'mi + 'ni2)ghs-i • Qi — W2g?2e3 • 92 for the gravitational acceleration g G M. Substituting these 
into (|2.10|) - ()2.1ip . the continuous equations of motion for the double spherical pendulum are given by 

qi=ujiy.qi q2=uj2'xq2, (3.19) 
(mi -I- 7712)^1/3x3 -m2lil2qiq2 

-m2lll2q2qi ™2^i/3x3 

which are more compact than existing equations written in terms of angles. Another nice property is that 
the same structure for the equations of motion is maintained for n > 2. Thus, it is easy to generalize these 
equations of motion to a triple, or more generally, a multiple-link spherical pendulum. 

We compare the computational properties of the discrete equations of motion given by (|3.10p - p.l3p 
with a 4(5)-th order variable step size Runge-Kutta method for (|3.19p - (|3.20p . We choose mi — m2 = 1kg, 
h = h ^ 9.81m. The initial conditions are 91, = [0.8660, 0, 0.5], g2o = [0, 0, 1], uji^ = [-0.4330, 0, 0.75], 
^2a = [0, 1, 0] rad/sec. The simulation time is 100 sec, and the step-size of the discrete equations of motion is 
h — 0.01. Figure [XT] shows the computed total energy and the configuration manifold errors. The variational 
integrator preserves the total energy and the structure of (§^)" well for this chaotic motion of the double 
spherical pendulum. The mean total energy variation is 2.1641 x 10~^Nm, and the mean unit length error is 
8.8893 X 10~^^. But, there is a notable increase of the computed total energy for the Runge-Kutta method. 





UJl 






UJ2 





7772^1^2(^2 • ^2)91 92 + {mi + m2)ghqie3 
m2/iZ2(t^i • wi)g2'Zi + 77725/29263 



(3.20) 
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t 

(b) Computed total energy 



20 40 60 80 

/ 

(c) Unit length error 



(a) Trajectory of pendulum 

Fig. 3.1. Numerical simulation of a double spherical pendulum (RK45: blue, dotted, VI: red, solid) 



where the mean variation of the total energy is 7.8586 x 10 Nm. The Runge-Kutta method also fails to 
preserve the structure of The mean unit length error is 6.2742 x 10 



-5 



Example 2 (n-body Problem on Sphere) An n-body problem on the two-sphere deals with the motion 
of n mass particles constrained to lie on a two-sphere, acting under a mutual potential. Let G M and 
qi G §^ be the mass and the position vector of the i-ih particle, respectively. The z, j-th element of the inertia 
matrix is My = rrii when i = j, and My = otherwise. In [S], the following expression for the potential is 
introduced as an analog of a gravitational potential. 



V{qi,. 



for a constant 7. Substituting these into (|2.7p . the continuous equations of motion for the n-body problem 
on a sphere are given by 



rrnqi 



-mi{qi ■ q.j)qi - qi x (qi x 7 



q] 



3 



1 (1 - (qi ■ qj) 



!)3/2 



(3.21) 



for i G {1, . . . , n}. 

A two-body problem on the two-sphere under this gravitational potential is studied in [6] by explicitly 
using unit length constraints. Here we study a three-body problem, n = 3. Since there are no coupling terms 
in the kinetic energy, we use the explicit form of the variational integrator. We compare the computational 
properties of the discrete equations of motion given by (|3.17p - (|3.18p with a 2-nd order fixed step size Runge- 
Kutta method for (j3.2ip . We choose mi = 7712 = = 1, and 7=1. The initial conditions arc qi^ = 
[0, -1, 0], g2o = [0, 0, 1], (73o = [-1, 0, 0], LOi„ = [0, 0, -1.1], lo2„ = [1, 0, 0], and = [0, 1, 0]. The 
simulation time is 10 sec. Figure [3^21 shows the computed total energy and the unit length errors for various 
step sizes. The total energy variations and the unit length errors for the variational integrator are smaller 
than those of the Runge-Kutta method for the same time step size by several orders of magnitude. For 
the variational integrator, the total energy error is reduced by almost 100 times from 1.1717 x lO^"' to 
1.1986 X IQ-s when the step size is reduced by 10 times from 10 ^ to 10 ^, which verifies the second order 
accuracy numerically. 
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(a) Trajectory of particles (b) Total energy error v.s. step size (c) Unit length error v.s. step size 

Fig. 3.2. Numerical simulation of a 3-body problem on sphere (RK2: blue, square, VI: red, circle) 



Example 3 (Interconnection of Spherical Pendula) We study the dynamics of n spherical pendula 
connected by linear springs. Each pendulum is a mass particle connected to a frictionless two degree-of- 
freedom pivot by a rigid massless link acting under a uniform gravitational potential. It is assumed that all 
of the pivot points lie on a horizontal plane, and some pairs of pendulua are connected by linear springs at 
the centers of links. 

Let the mass and the length of the i-th pendulum be mi^U G R, respectively. The vector qi G 
represents the direction from the i-ih. pivot to the i-th mass. The inertia matrix is given by Mij = mill 
when i = j, and = otherwise. Let S be a set defined such that G S if the i-th pendulum and 

the j-th pendulum are connected. For a connected pair G S, define Kij G M and r^j G M.^ as the 

corresponding spring constant and the vector from the i-th pivot to the j-th pivot, respectively. The bases 
for the inertial frame are chosen such that the direction along gravity is denoted by 63 = [0, 0, 1] G M"^, and 
the horizontal plane is spanned by ei = [0, 0, 1], 62 = [0, 1, 0] G K."^. The potential energy is given by 




Substituting these into (|2.9|) - (|2.10p . the continuous equations of motion for the interconnection of spherical 
pendula are given by 

dV 

rriilfLO, = -q,x— (3.22) 
oqi 

qi = LOi X qi (3.23) 

for i G {1, . . . , n}. 

We compare the computational properties of the discrete equations of motion given by (|3.17[) - (|3.18[) 
with a 2-nd order fixed step size explicit Runge-Kutta method for p.22p - p.23p . and the same Runge-Kutta 
method with reprojection; at each time step, the vectors qi^ are projected onto §^ by using normalization. 

We choose four interconnected pendula, n = 4, and we assume each pendulum has the same mass and 
length; nii = 0.1kg, li — 0.1m. The pendula are connected as S = {(1, 2), (2, 3), (3, 4), (4, 1)}, and the 
corresponding spring constants and the relative vector between pivots are given by K12 = 10, K12 — 20, 
K12 = 30, K12 — 40N/m, ri2 = —^34 — kei, and r23 = —rn = —lie2- The initial conditions are chosen as 
qi„ = <72q = <74o = 63, q3„ ^ [0.4698,0.1710,0.8660], llIi„ = [-10,4,0], and uj2a = ^Sq = = Orad/sec 
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(a) Motion of pendula (b) Computed total energy (c) Unit length error 

Fig. 3.3. Numerical simulation of a system of 4 spherical pendula (RK2: blue, dotted, RK2 with projection: black, dashed, 
VI: red, solid) 

Figure 13.31 shows the computed total energy and the unit length errors. The variational integrator 
preserves the total energy and the structure of (S^)" well. The mean total energy variation is 3.6171 x 
10~^Nm, and the mean unit length error is 4.2712 x 10^^^. For both Runge-Kutta methods, there is a 
notable increase of the computed total energy. It is interesting to see that the reprojection approach makes 
the total energy error worse, even though it preserves the structure of accurately. This shows that a 

standard reprojection method can corrupt numerical trajectories [51 113j. 



Example 4 (Pure Bending of Elastic Rod) We study the dynamics of {n + 1) rigid rod elements that 
are serially connected by rotational springs, where the 'zeroth' rod is assumed to be fixed to a wall. Thus, 
the configuration space is (§^)". This can be considered as a simplified dynamics model for pure non-planar 
bending of a thin elastic rod that is clamped at one end and free at the other end. Notably, this approach 
is geometrically exact, and preserves the length of the elastic rod in the presence of large displacements. 

The mass and the length of the i-th rod element are denoted by mi,li e R, respectively. The inertia 
matrix is given by 

n n ^ 

k—i-\-l A:— max{2,j } 

for i,j G {1, . . . n} and i ^ j. The potential energy is composed of gravitational terms and elastic bending 
terms given by 

n i—1 ^ -j^ 

V{qi,.. . ,g„) = -^mig{^ljqj + -kqi) ■ 63 + -k^(1 - qi-i ■ qif, 
i=l j=l 

where a constant vector go £ denotes the direction of the zeroth rod element fixed to a wall, and G M 
denotes spring constants. The bases for the inertial frame are chosen such that the gravity direction is 
denoted by 63 = [0,0,1] £ R'^, and the horizontal plane is spanned by ei = [0,0, 1], 62 = [0,1,0] G R'^. 
Suppose that the total mass and length of rod are given by m, I, and each rod element has the same mass 
and length, i.e. nii = li = for i e {0, ...,n}. Substituting these into (|2.8p . the continuous 
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equations of motion for the pure bending of an elastic rod are given by 

~ («+l)^ ™ (92-92)92+92 957 

(3.24) 

We compare the computational properties of the discrete equations of motion given by (|3.10p - (|3.13|) 
with a 4(5)-th order variable step size Runge-Kutta method for (|3.24p . We choose 10 rod elements, n = 10, 
and the total mass and the total length are m = 55 g, I = 1.1, m. The spring constants are chosen as 
Ki = 1000 Nm. Initially, the rod is aligned horizontally; = ei for all i G l,...rt. The initial angular 
velocity for each rod element is zero except 0^5^ — [0, 0, 10] rad/sec. This represents the dynamics of the rod 
after an initial impact. The simulation time is 3 sec, and the step-size of the discrete equations of motion is 
h = 0.0001. 

Figure 13.41 shows the computed total energy and the unit length errors. The variational integrator 
preserves the total energy and the structure of (S^)". The mean total energy variation is 1.4310 x 10^^ Nm, 
and the mean unit length error is 2.9747 x 10^^**. There is a notable dissipation of the computed total 
energy for the Runge-Kutta method, where the mean variation of the total energy is 3.5244 x 10~*Nm. The 
Runge-Kutta method also fail to preserve the structure of (§^)". The mean unit length error is 1.8725 x 10~^. 



(n+l)^ 
n-1 



3x3 



n-1 



-2(;7^TO^'9n9n 



2ii^ml 9191 

n— 5/3 72 T 



2(n+l)''"^^^9n9ri 



91 




92 




.9n. 





Example 5 (Spatial Array of Magnetic Dipoles) We study dynamics of n magnetic dipoles uniformly 
distributed on a plane. Each magnetic dipole is modeled as a spherical compass; a thin rod magnet supported 
by a frictionless, two degree-of-freedom pivot acting under their mutual magnetic field. This can be considered 
as a simplified model for the dynamics of micromagnetic particles [3]. 

The mass and the length of the i-th magnet are denoted by mi,li G M, respectively. The magnetic 
dipole moment of the i-th magnet is denoted by ViQi, where j/^ G M is the constant magnitude of the 
magnetic moment measured in ampere square-meters, and qi G is the direction of the north pole from 
the pivot point. Thus, the configuration space is (S^)". The inertia matrix is given by Mij = t-^toJ| when 
i — j, and Mij — otherwise. Let G M"^ be the vector from the i-pivot point to the j-th pivot point. The 
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Fig. 3.5. Numerical simulation of an array of magnetic dipoles (RK45: blue, dotted, VI: red, solid) 



mutual potential energy of the array of magnetic dipoles are given by 



1 " 

V{qi,...,qn) ^ l^Yl 



(it -Qj) 



where /J, = 47rx 10 ^N-A ^ is the permeability constant. Substituting these into (|2.9p - (|2.10p . the continuous 
equations of motion for the spatial array of magnetic dipoles are given by 



1 " 
—rmlfCo, = -q, X ^ 



1j 



qi = UJi X qi 



(3.25) 
(3.26) 



for i e {1, . . . , n}. 

We compare the computational properties of the discrete equations of motion given by (|3.17p - (|3.18p with 
a 4(5)-th order variable step size Runge-Kutta method for (|3.25p - (|3.26p . We choose 16 magnetic dipoles, 
n = 16, and we assume each magnetic dipole has the same mass, length, and magnitude of magnetic moment; 
rrii = 0.05 kg, U = 0.02 m, Vi = 0.1 A • m^. The magnetic dipoles are located at vertices of a 4 x 4 square grid 
in which the edge of a unit square has the length of 1.2Zi. The initial conditions are chosen as qig = [1, 0, 0], 
w,, = [0, 0, 0] for alH G {1, . . . , 16} except qie^ = [0.3536, 0.3536, -0.8660] and wi„ = [0, 0.5, 0] rad/sec. 

Figure 13.51 shows the computed total energy and the unit length errors. The variational integrator 
preserves the total energy and the structure of (S^)" well. The mean total energy variation is 8.5403 x 
lO"!" Nm, and the mean unit length error is 1.6140 x 10"^''. There is a notable dissipation of the computed 
total energy for the Runge-Kutta method, where the mean variation of the total energy is 2.9989 x 10~''Nm. 
The Runge-Kutta method also fail to preserve the structure of (§^)". The mean unit length error is 1.7594 x 
10"-*. 



Example 6 (Molecular Dynamics on a Sphere) We study molecular dynamics on §-^. Each molecule 
is modeled as a particle moving on EP. Molecules are subject to two distinct forces: an attractive force at 
long range and a repulsive force at short range. Let rrii G M and qi £ E>^ be the mass and the position 
vector of the i-th molecule, respectively. The «, j-th element of the inertia matrix is Mij = mi when i — j, 
and Mij — otherwise. The Lennard- Jones potential is a simple mathematical model that represents the 
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(a) i = 



(a) Initial trajectories (b) Computed total energy 

Fig. 3.6. Numerical simulation of molecular dynamics on a sphere 




(b) t = 0.25 (c) t = 0.5 (d) t = 0.75 

Fig. 3.7. Kinetic energy distributions over time 



(e) t = 5 



\( 


r- 


( a 













behavior of molecules [12] 

1 " 

y((Zi,...,g„) = 2 II 

where the first term models repulsion between molecules at short distance according to the Pauli principle, 
and the second term models attraction at long distance generated by van der Walls forces. The constant e 
and a are molecular constants; e is proportional to the strength of the mutual potential, and cr characterize 
inter- molecular force. Substituting these into (|2.7p . the continuous equations of motion for the molecular 
dynamics on a sphere are given by 



miiii = ~m,{q, ■ qi)q, - q, x [qi x y^4e ^' 



11% -qjW 



via 



12 



(3.27) 



for i e {1, . . . , n}. 

We choose 642 molecules, n — 642, and we assume each molecule has the same mass, nii = 1. Initially, 
molecules are uniformly distributed on a sphere. The strength of the potential is chosen as e = 0.01, and 
the constant a is chosen such that the inter-molecular force between neighboring molecules is close to zero. 
The initial velocities are modeled as two vortices separated by 30 degrees. The simulation time is 5 sec, and 
the step size is h = 0.005. 

Trajectories of molecules and the computed total energy is shown at Figure [3T6l The mean deviation of 
the total energy is 1.8893 x 10"'^, and the mean unit length error is 5.2623 x 10~^^. In molecular dynamics 
simulations, macroscopic quantities such as temperature and pressure are more useful than trajectories 
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of molecules. Figure 13.71 shows the change of kinetic energy distributions over time, which measures the 
temperature [l]; the sphere is discretized by an icosahedron with 5120 triangular faces, and the color of 
a face is determined by the average kinetic energy for molecules that lie within the face and within its 
neighboring faces. The local kinetic energy is represented by color shading of blue, green, yellow, and red 
colors in ascending order. 

4. Conclusions. Euler-Lagrange equations and variational integrators are developed for Lagrangian 
mechanical systems evolving on (S^)" where the Lagrangian is written in a particular form given by (|2.2p . 
The structure of §^ is carefully considered to obtain global equations of motion on (§^)" without local 
parameterizations or explicit constraints. 

In the continuous time setting, this provides a remarkably compact form of the equations of motion 
compared to the popular angular description. For example, it is not practical to study a triple spherical 
pendulum by using angles due to the complexity of the trigonometric expressions involved. The global 
Euler-Lagrange equations on (S^)" maintain the same compact structure for arbitrary n. In particular, it is 
possible to use them as a finite element model for a continuum problem as shown in Example (4] They are 
also useful for the theoretical study of global dynamic characteristics. 

The variational integrators on (§^)" preserve the geometric properties of the dynamics as well as the 
structure of the configuration manifold concurrently. They are symplectic, momentum preserving, and they 
exhibit good energy behavior for exponentially long time as they are derived from discrete Hamilton's prin- 
ciple. Using the characteristics of the homogeneous manifold (§^)", the discrete update map is represented 
by a group action of S0(3) to obtain compactly represented, and possibly explicit, numerical integrators. In 
particular, variational integrators on (§^)" completely avoid the singularities and complexity introduced by 
local parameterizations and explicit constraints. 
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